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I .    INTRODUCTION 

A.   BACKGROUND 

Autonomous  vehicles  suitable  for  use  in  modern 
applications  require  high  maneuverability,  quick  response,  and 
robust  performance  characteristics  [Ref.  7].  In  order  to 
maintain  accurate  track  keeping  characteristics,  a  suitable 
combination  of  path  planning,  navigation,  guidance,  and 
control  design  is  needed  [Ref.  10]  .  Sufficient  information  is 
obtained  from  charted  obstacles  and  the  environment  for  smooth 
paths  to  be  generated  for  the  vehicle  to  follow  [Ref.  8]. 
Although  it  is  possible  to  design  a  nonlinear  controller  which 
incorporates  and  utilizes  information  on  the  nonlinear  dynamic 
properties  of  the  vehicle,  as  well  as  the  geometric 
nonlinearities  of  the  desired  track,  the  resulting  scheme 
tends  to  be  very  complex  and  time  consuming  [Ref.  3].  The 
alternative  is  to  use  separated  navigation,  guidance,  and 
autopilot  functions.  A  certain  level  of  feedback  is  provided 
at  the  navigation  level  through  the  use  of  sonars  in  order  to 
replan  a  path  due  to  uncharted  obstacles  or  changes  in  mission 
objectives.  This  operation  is  event-driven  and  occurs 
asynchronously,  and  in  many  cases  it  does  not  affect  stability 
and  performance  of  the  overall  navigation,  guidance,  and 
control   scheme.     Based   on   the   provided   navigational 


information,  the  guidance  law  generates  heading  commands, 
which  are  executed  by  the  control  law  by  appropriate  use  of 
vehicle  effectors,  such  as  control  surfaces  and  cross  body 
thrusters.  However,  the  time  required  to  process  sonar  data 
and/or  inertial  position  information  may  be  significant  and 
cannot  be  neglected  [Ref .  9] .  In  addition,  the  guidance  and 
control  laws  must  be  as  fast  as  possible  in  order  to  ensure 
accurate  path  keeping  characteristics.  The  navigational 
position  information  time  lags,  as  well  as  the  dynamic  lags 
due  to  the  vehicle  inertia,  however,  set  a  limit  on  the 
vehicle  reaction  time.  Therefore,  stability  of  the  combined 
scheme  becomes  an  issue  which  requires  analysis. 

B.   OBJECTIVE 

Previous  studies  have  established  boundaries  for  guidance 
and  control  laws  in  the  horizontal  [Ref.  10]  and  vertical 
plane  [Ref.  12]  maneuvering  along  straight  line  paths,  as  well 
as  curved  segments  [Ref.  11].  The  most  important  assumption 
in  those  results  was  the  existence  of  instantaneous  positional 
information  updates  when  needed.  In  this  work  we  relax  the 
latter  assumption.  Stability  boundaries  are  computed 
parametrized  by  the  amount  of  positional  information  time  lag. 
Results  are  presented  based  on  eigenvalue  and  frequency 
response  techniques.  This  thesis  builds  on  the  linear 
stability  analysis  performed  in  [Ref.  13]  which  can  predict 
the  shape  of  the  stability  boundaries.   In  this  work,  vehicle 


motions  in  the  vicinity  of  the  corresponding  boundary  are 
assessed  using  nonlinear  bifurcation  theory  [Ref .  5] .  It  is 
shown  that  the  loss  of  stability  occurs  as  generic  Hopf 
bifurcations,  where  upon  loss  of  stability  of  straight  line 
motion  a  family  of  periodic  solutions  appears.  Nonlinear 
expansions  utilizing  center  manifold  approximations  and 
integral  averaging  techniques,  reveal  that  these  Hopf 
bifurcations  can  occur  either  in  supercritical  or  subcritical 
forms.  In  the  supercritical  case,  a  stable  family  of  limit 
cycles  is  generated  immediately  after  the  nominal  motion 
becomes  unstable.  In  the  subcritical  case,  however,  the 
resulting  limit  cycles  are  initially  unstable  and  they  are 
generated  even  before  the  nominal  motion  loses  its  stability. 
In  the  latter  case,  the  domain  of  convergence  of  straight  line 
motion  becomes  increasingly  smaller  as  the  stability  boundary 
is  reached.  As  a  result,  the  methods  developed  in  this  work 
characterize  the  degree  of  confidence  of  the  computed 
stability  boundaries  by  examining  stability  under  finite 
external  disturbances.  A  first  order  memory  scheme,  which  can 
be  easily  implemented  in  real  time,  is  suggested  in  order  to 
expand  the  region  of  stability  of  straight  line  motion.  All 
computations  are  performed  for  the  NPS  autonomous  underwater 
vehicle  for  which  a  set  of  geometric  properties  and  slow 
motion  hydrodynamic  derivatives  is  available  [Ref.  1]  .  Unless 
otherwise  mentioned,  all  results  are  presented  in  standard 
nondimensional   form.    Equation  development  presented  in 


[Ref .13]  has  been  used  in  chapter  II,  and  with  modifications 
in  chapter  IV. 

C.   THESIS  OUTLINE 

Chapter  II  presents  the  formulation  of  vehicle  equations 
of  motion  and  the  criterion  for  determining  the  region  of 
guidance/control  stability  in  the  presence  of  a  position 
information  time  lag. 

Chapter  III  examines  the  vehicle's  control  stability 
through  the  use  of  a  Hopf  bifurcation  analysis. 

Chapter  IV  examines  the  effect  on  the  region  of  stability 
determination  by  the  use  of  a  memory  model  which  incorporates 
the  two  previous  vehicle  positions. 


II.   PROBLEM  FORMULATION 

A.   EQUATIONS  OF  MOTION 

The  mathematical  model  which  describes  vehicle  motion  in 
the  horizontal  plane  consists  of  nonlinear  sway  and  yaw 
equations  of  motion.  A  moving  coordinate  frame  which  is  fixed 
at  the  vehicle's  geometric  center  gives  Newtonian  equations 
of  motion  which  are 

m(v+ur+xGf-yGr2)  =Y  (2.1) 

Izr  +mxG{  v+ur)  -myGvr  =N  (2.2) 

where  u  is  constant  forward  speed;  v  and  r  are  the  relative 
sway  and  yaw  velocities  of  the  moving  vehicle  relative  to  the 
water;  m  is  the  mass  of  the  moving  vehicle;  Xq,  yG  are  the 
respective  lateral  and  longitudinal  positions  of  the  center  of 
gravity;  and  Y,  N  represent  the  total  excitation  sway  force 
and  yaw  moment,  respectively.  The  vehicle's  added  mass, 
damping,  and  viscous  drag  due  to  motion  through  the  water  are 
represented  by  quadratic  drag  terms  and  memoryless  polynomials 
in  v  and  r.  Y  and  N  can  be  represented  as  the  sum  of  these 
terms  by 


Y=£liYir  +  -£-l*{Yvv+YIur)  +-£-l2Yvuv 


N=£l*Ni±+£l*(Ni,v+Nziir)  +-^-l2Nvuv 


where  1  is  the  vehicle  length,  Ya,  Na  represent  partial 
derivatives  of  Y  and  N  with  respect  to  a,  CDy  is  the  drag 
coefficient,  and  6    is  the  rudder  angle. 

At  relatively  high  speeds  and  low  angles  of  attack  the 
steering  response  is  predominantly  linear.  In  low  speed 
maneuvering  or  hovering  conditions  the  cross  flow  integral 
drag  term  becomes  significant. 

Vehicle  yaw  rate  is  expressed  by 

i|r  =  r  (2.3) 

and  inertial  position  rate  is  expressed  by 

y  =  usinijj  +  vcosijf  (2.4) 

where  \p   is  the  heading  angle  and  shown  in  Figure  1. 


Figure  1.   Vehicle  Geometry  and  Definition  of  Symbols. 


Taking  equations  (2.1),  (2.2),  and  (2.3)  as  a  set  of  three 
coupled  linear  differential  equations,  and  a  linearized  form 
of  equation  (2.4),  the  final  set  of  equations  of  motion  for 
steering  control  are 

i|r=r  (2.5a) 

v=a11uv+a12ur+b1u2b  (2.5b) 


r  =  a21uv+a22ur+b2u2b  {2.5c] 


y=ui|r  +  v  (2.5d) 

at  any  nominal  forward  speed. 

B.   GUIDANCE  AND  CONTROL 
1.   Guidance  Scheme 

An  orientation  based  command  scheme  is  used  where  the 
commanded  vehicle  heading  angle  \pc  is  a  function  of  the  line 
of  sight  angle  o  between  the  actual  vehicle  position  and  a 
reference  position  on  the  nominal  path  located  at  a  constant 
lookahead  distance  d.  Schematically,  this  is  presented  in 
Figure  1 . 

The  line  of  sight  angle  a  is  defined  by 


tana  =  --^  (2.6) 


The  autopilot  will  be  called  upon  to  orient  the 
vehicle  to  the  commanded  heading  angle  \pc,  which,  in  pursuit 
guidance,  will  equal  the  line  of  sight  angle.  This  defines  \pc 
as 


i|ic=-arctan^  (2.7) 


For  relatively  small  angles 


♦.""§ 


2.   Controller  Design 

The  steering  control   governing  equations   can  be 
represented  in  the  form 


x=Ax+bb  (2.9) 

where  the  state  vector  equation  is 

x=[x\>,v,r]T  (2.10) 

Linear  full  state  feedback  is  introduced  in  the  form 


6=A:1(^-i|rc)  +k2v+k3r  (2.11) 

The  gains  k1#  k2 ,  and  k3  are  computed  by  pole  placement 
to  give  the  desired  dynamics  to  the  closed  loop  system.  The 
vehicle's  longitudinal  axis  is  pointed  toward  the  desired 
heading  by  the  difference  (\J/-\pc)    in  the  control  law. 

With  a  dimensionless  time  constant  tc  the  general  form 
of  the  characteristic  equation  is 

X3+a1X2  +  a2A.*a3  =  0  (2.12) 


where 


±,a2  =  ^,a,  =  ± 
t'2      t2       3   t3 


The  controller  gains  are  computed  from 


(b2a11-b1a21)  u 


(bxa22-b2a12)  u2k2  +  (b2axx-bxa21)  u3k2 
=  a2+b2u2kx-  (axxa22-ax2a2X)  u2 


2.13a: 


(2.13b) 


bxu2k2+b2u2k2  =  -ax  -  (axl+a22)  u  (2  .  13c) 
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These  gain  values  are  continuously  updated  due  to  changing 
nominal  forward  speeds.  Figure  2  shows  the  three  controller 
gains  at  unit  forward  speed  versus  the  system  time  constant. 

C.   TIME  LAG 

All  of  the  required  state  information  for  vehicle  control 
is  available  at  sufficient  rates  with  the  possible  exception 
of  the  y  position  information.  Due  to  time  requirements  for 
reduction  of  sonar  returns  and  inertial  navigation 
information,  there  may  be  a  time  lag  in  the  y  position 
information. 

This  time  lag  T  (sec)  is  introduced  to  the  guidance 
control  law  and  the  commanded  rudder  angle  \pc   becomes 


i|rc=-arctany(t  r)  (2  .  14) 


For  small  angles 


tc=-zi^l  (2.i5: 


and  the  resulting  linearized  control  law  becomes 

b=k1i\>+k1y{t,T)  +k2v+k3r  (2.16) 
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Figure  2.   Controller  Gains  versus  tc  at  unit  forward  speed, 
(k^solid,  k2=dashed,  k3=dotted) 
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The  linearized  equations  of  motion  become 


2.17a: 


v  =  b1u2k1i\j+  (a11u+jb1u2ic2)  v 

k,  (2.17b) 

+  ia^u+b^u2^)  r+Jb1u2-iy(  t-T) 


z=b2u2k1^  +  (a21u+b2u2k2)  v 

, K      .         .  (2.17c) 


+  (a22u+Jb2u2Jc3)  r+b2u2-±y(t-T) 

y=Wi\i  +  v  (2.17d) 

D.   STABILITY  ANALYSIS 

Control  law  stability  is  determined  by  the  eigenvalues  of 
the  matrix  A  from  the  linear  system 

x  =  Ax  (2.18) 

with  the  state  vector 

x=[y\>,v,r,y]T  (2.19) 
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where  the  state  equations  have  been  linearized  around  straight 
line  motion.  The  characteristic  equation  of  (2.18)  is  a 
quartic  of  the  form 

V+BV+C\2+DX+E=0  (2.20) 

where  the  coefficients  B,  C,  D,  and  E  are  functions  of  vehicle 
properties,  controller  gains,  and  lookahead  distance. 

The  characteristic  equation  will  give  a  pair  of  complex 
conjugate  roots  which  cross  the  imaginary  axis  under  the 
following  conditions. 

BCD-B2E-D2=0  (2.21a) 


—  >0  (2.21b) 

B 


Minimum   lookahead   distance   dcrit   for   stability   can   be 
determined  by  translating  these  conditions  to 

a1d2+a2d+a3=0  (2.22a) 


jb,  a5,-jb,a1  ~-b0 

d>   — i-22 212 2  (2.22b) 


where 
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a1=a1a2-a2  (2.23a) 

(a1a2-2oc3)  (b1a22-b2a12-b2) 

(2.23b) 


(b2axl-bra21)  u 


_  -  {bxa22-b2a12-b2)  [b.a^  (b±a22-b2a12-b2)  u]  a3 

<33 (2.23c 

{b2a11-b1a21)2u 


Analysis  of  the  system  when  operating  with  an  information 
time  lag  requires  approximation  by  either  a  first  order  Taylor 
expansion 

y(t-T)  =y-Ty  (2.24) 

a  second  order  Taylor  expansion 

y(t-T)  =y-Ty+—y  (2.25) 


a   third  order  Taylor  expansion 


rn2  rp3 


y(t-T)  =y-Ty+  —  y-  —  y  (2.26) 

2  6 


or  a  frequency  response  analysis  using  the  Nyquist  criterion. 
Figure  3  shows  the  region  of  stability  given  by  each  of  these 
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methods  with  a  one  second  (dimentionless)  time  lag.  It  can  be 
seen  that  the  agreement  is,  in  general,  very  good  among  the 
three  Taylor  series  approximations,  especially  as  the  time 
constant  tc  is  increased.  The  third  order  approximation 
coincides  with  the  exact  solution  from  the  frequency  response 
(Nyquist  criterion)  method.  The  first  order  method  requires 
the  highest  value  of  d  for  stability  at  a  given  tc,  and 
therefore  is  the  most  conservative  method  for  design  use.  For 
this  reason  the  first  order  method  is  used  in  the  next  chapter 
to  study  the  dynamics  of  the  system  as  the  critical  stability 
boundary  is  crossed. 
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I    1.2  - 


0.8 


2.H 


Figure  3.  Taylor  expansion  and  Nyquist  stability  analysis 
with  a  one  second  time  lag.  (Taylor  first 
order=solid,  Taylor  second  order=dashed,  Taylor 
third  order  and  Nyquist=dotted) 
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III.   HOPF  BIFURCATION 

A.   INTRODUCTION 

It  can  be  numerically  confirmed  that  in  all  cases  of 
stability  loss  of  the  previous  chapter,  one  pair  of  complex 
conjugate  eigenvalues  of  the  corresponding  eigenvalue  problem 
crosses  transversely  the  imaginary  axis.  A  situation  like 
this  in  which  a  certain  parameter  is  varied  such  that  the  real 
part  of  one  pair  of  complex  conjugate  eigenvalues  of  the 
linearized  system  matrix  crosses  zero,  will  result  in  the 
system  leaving  its  steady  state  in  an  oscillatory  manner. 
This  loss  of  stability  is  called  Hopf  bifurcation  and 
generically  occurs  in  one  of  two  ways,  supercritical  or 
subcritical.  In  the  supercritical  case,  stable  limit  cycles 
are  generated  after  the  nominal  straight  line  motion  loses  its 
stability.  The  amplitudes  of  these  limit  cycles  are 
continuously  increasing  as  the  parameter  distance  from  its 
critical  value  is  increased.  For  small  values  of  this 
criticality  distance  the  resulting  limit  cycle  is  of  small 
amplitude  and  differs  little  from  the  initial  nominal  state. 
In  the  subcritical  case,  however,  unstable  limit  cycles  are 
generated  before  the  nominal  state  loses  its  stability. 
Therefore,  depending  on  the  initial  conditions  it  is  possible 
to  diverge  away  from  the  nominal  straight  line  path  and 
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converge  toward  a  different  steady  state  even  before  the 
nominal  motion  loses  its  stability.  In  many  cases  this  new 
steady  state  is  another  limit  cycle  of  considerably  higher 
amplitude.  This  means  that  in  the  subcritical  Hopf 
bifurcation  case  the  domain  of  attraction  of  the  nominal  state 
is  decreasing  and  in  fact  it  shrinks  to  zero  as  the  critical 
point  is  approached.  Random  external  disturbances  of 
sufficient  magnitude  can  throw  the  vehicle  off  to  an 
oscillatory  steady  state  even  though  the  nominal  state  may 
still  remain  stable.  After  the  nominal  state  becomes 
unstable,  a  discontinuous  increase  in  the  magnitude  of  motions 
is  observed  as  there  exists  no  simple  stable  nearby  attractors 
for  the  vehicle  to  converge  to.  Distinction  between  these  two 
qualitatively  different  types  of  bifurcation  is,  therefore, 
essential  in  design.  The  computational  procedure  requires 
higher  order  approximations  in  the  equations  of  motion  and  is 
the  subject  of  the  next  section. 

B .   COMPUTATIONS 

The  vehicle  equations  of  motion  are  written  in  the  form 

i|i  =  r  (3.1a) 


v=a11uv+a12ur+b1u2b  (3.1b) 
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f  =  a21uv+a22ur+b2u2S  (3.1c) 


y=usini|/  +  vcosi|r  (3. Id) 


In  order  to  capture  the  effect  of  rudder  saturation  we 
write  the  rudder  angle  6  in  the  form 

6=6sat  tanM-jr-^-)  (3.2) 

"sat 


where  <5sat  is  the  saturation  limit  of  <5 ,  typically  set  at  +0.4 
radians.  Equation  (3.2)  is  used  instead  of  a  hard  saturation 
function  because  of  its  smoothness,  which  is  important  for  the 
computation  in  this  section.  60  is  the  linear  control  law  of 
the  form 


b^^iq-qj+^v+k^r  (3.3 


where 


tamttc=-y(t:r)  —^P-  (3.4) 


using  the  first  order  approximation  for  y(t-T). 

The  state  equations  (3.1)  through  (3.4)  are  written  in  the 
compact  form 
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x=f(x)  (3.5) 

where  f  (x)  is  a  nonlinear  function  of  the  state  variables 
vector 

x=  [i|i,  v,i,y]  (3.6) 

Expanding  f  (x)  in  Taylor  series,  we  can  write  (3.5)  in  the 
form 


x  =  Ax+g(x)  (3.7) 

where  A  is  the  linear  system  matrix  and  g (x)  contains  all  the 
leading  nonlinear  terms  of  f (x) .  Due  to  the  port/starboard 
symmetry  of  our  problem,  g(x)  contains  only  third  order  terms. 
Defining  T  as  the  matrix  of  eigenvalues  of  A  and  applying  the 
transformation 

x=Tz  (3.8) 

the  state  equation  is  transformed  into  its  canonical  form 

z  =  T-lATz  +  T~1g{Tz)  (3.9) 

At  the  Hopf  bifurcation  point 
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where  p  and  q  are  negative  and  co0  is  positive, 
system  of  coordinates 


(3.10) 


In  the  new 


z  =  T~xx 


(3.11) 


the  dynamics  of  the  system  are  generated  by  a  reduced  two 
dimensional  system  z1  and  z2 ,  since  the  coordinates  z3  and  z4 
correspond  to  the  eigenvalues  p  and  q  and  are  asymptotically 
stable.  These  stable  variables  z3 ,  z4  can  be  expressed  as  a 
function  of  the  critical  variables  zlt  z2;  these  functional 
expressions  are  at  least  of  third  order.  Therefore,  the 
stable  coordinates  z3,  z4  do  not  influence  the  third  order 
expansion  of  (3.7) .  Using  the  above  expression,  we  can  write 
the  two  dimensional  system  in  z±,    z2  as 


*1  =  -UoZ^I^zl  ^T^zlz-^I^Z^l^I^zl 


(3.12a; 


2T2  =W0Zl+r21Zl+r22Zl2*2+r23*lZ22+-r24Z2 


(3.12b) 


where  the  coefficients  r-  are  complicated  expressions  that 
are  derived  from  (3.9) 
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Equations  (3.12)  are  valid  at  exactly  the  Hopf 
bifurcation  point.  For  values  of  the  parameter  d  close  to  the 
bifurcation  point,  they  are  written  as 

zx  =  a'cz1-  (Wo+co'e)  z2+Fx  (z1,  z2)  (3  .  13a) 


z2  =  (cOq+co^)  zl  +  alzz2+F2  (z1,  z2)  (3  .  13b) 

where  ol  ',  to  '  are  the  derivatives  of  the  real  and  imaginary 
parts  of  the  critical  pair  of  eigenvalues  evaluated  at  the 
bifurcation  point;  e  is  the  difference  of  the  parameter  d  from 
its  critical  value;  and  the  nonlinear  function  F1#  F2  remain 
the  same  as  in  (3.12), 


F1(zllz2)  =r11zl+r12z?z2+r12z1z2+r14z2 


If  we  introduce  polar  coordinates  in  the  form 


equations  (3.13)  become 


(3.14a; 


F2(z1,z2)  =r21zl+r22z?z2+r23z1z;+r24Z2  (3.14b! 


(3.15a! 


(3.15b: 
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R  =  a'zR+F1(R,Q)  cos6+F2(i?,8)  sin0  (3.16a) 

i?6  =  ((oo+a)/e)i?+F2(i?,6)cose-F1(i?,e)sin0  (3.16b) 

Equation    (3.16a)    yields 

R  =  a/eR+P(Q)R3  (3.17) 

where  P  {6+2tt)  =P  ( 6)   .   If  (3.17)  is  averaged  over  one  cycle  in 
6,    we  get  an  equation  with  constant  coefficients 

R  =  a'eR+KR3  (3.18) 


where 


K=^-  f2*P(Q)dQ  (3.19) 

2n  Jo 


Carrying  out  the  indicated  integration  in  (3.19)  we  obtain 
K=±(2r11  +  r13+r22  +  3r24)  (3.20 


Existence  and  stability  of  limit  cycles  can  be  determined 
by  analyzing  the  equilibrium  points  of  the  averaged  equation 
(3.18),  which  correspond  to  periodic  solutions  in  z1#  z2  as 
can  be  seen  from  (3.15)  .  If  we  examine  equation  (3.18)  we  can 
see  that : 
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1)  If  a   >0,  then: 

(a)  If  K>0  then  unstable  limit  cycles  coexist  with  the 
stable  equilibrium  for  £<0. 

(b)  If  K<0  then  stable  limit  cycles  coexist  with  the 
unstable  equilibrium  for  £>0. 

2)  If  a  '<0,  then: 

(a)  If  K>0  then  unstable  limit  cycles  coexist  with  the 
stable  equilibrium  for  £>0. 

(b)  If  K<0  then  stable  limit  cycles  coexist  with  the 
unstable  equilibrium  for  £<0. 


Therefore,  we  can  see  that  computation  of  the  above 
nonlinear  coefficient  K  can  distinguish  between  the  two 
different  types  of  Hopf  bifurcation: 

•  Supercritical  if  K<0 

•  Subcritical  if  K>0 
[Ref .  10] . 

C .   RESULTS 

Values  of  the  coefficient  K  have  been  calculated  for 
several  different  values  of  a  y  position  information  time  lag, 
over  a  span  of  system  time  constants,  using  the  Fortran 
program  presented  in  appendix  A.  This  program  has  also  been 
used  to  determine  the  vehicle  period  of  oscillation  in  the 
same  conditions.  Figure  4  shows  the  behavior  of  the  K 
coefficient  at  three  different  time  lags,  and  Figure  5  shows 
periods  of  oscillation  at  the  same  three  time  lags.   It  can  be 
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I.M 


Figure   4.      Coefficient  K  versus  tc  at  three  values  of  time  lag 
T.       (0.5    sec=solid,    1.0    sec=dashed#    1.5    sec=dotted) 
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1.13 


Figure  5.  Vehicle  oscillation  period  versus  tc  for  three 
values  of  time  lag  T.  (0.5  sec=solid,  1.0 
sec=dashedf  1.5  sec=dotted) 
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seen  that  supercritical  bifurcations  are  ensured  for 
sufficiently  high  values  of  the  control  time  constant  tc.  For 
tc  less  than  a  certain  critical  value,  the  corresponding  Hopf 
bifurcations  are  subcritical.  This  situation  should  be 
avoided  in  practice. 

D.   SIMULATIONS 

The  vehicle  path  has  been  simulated  by  using  the  Euler 
integration  method  coded  in  a  Fortran  program  presented  in 
appendix  B.  The  vehicle  control  law  (2.17)  as  well  as  the 
control  gains  (2.13)  are  used.  The  program  has  been  run  with 
input  values  for  system  time  constant,  y  position  information 
time  lag,  and  lookahead  distance.  The  vehicle  was  given  an 
initial  nominal  y  offset  of  half  a  shiplength  and  a  nominal 
forward  speed.  The  resulting  plots  of  y  position  against  time 
show  the  vehicle's  oscillatory  path  either  converge  to  the 
commanded  path  or  diverge. 

The  results  of  these  simulation  runs  were  compared  to  the 
results  of  the  Hopf  bifurcation  computations  in  two  ways.  The 
period  of  oscillation  observed  in  the  vehicle  path  was 
compared  to  the  period  predicted  by  the  Hopf  program. 
Additionally,  the  vehicle  stability,  determined  by  convergence 
to  the  commanded  path,  was  compared  to  the  K  coefficient  sign 
prediction  of  stability  which  was  determined  by  the  Hopf 
computations.  Figure  6  shows  a  vehicle  path  in  stable 
conditions.    Figure  7  shows  a  vehicle  path  in  unstable 
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conditions.  These  simulations  were  chosen  in  the 
supercritical  bifurcation  case  and  it  can  be  seen  that  the 
period  of  oscillation  observed  from  Figures  7  agrees  very  well 
with  the  theoretical  results  of  Figure  5. 
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J  00 


Figure  6.  Vehicle  path  in  stable  conditions  where  time 
constant  tc  =  1.0  sec,  time  lag  T  ■  1.0  sec,  and 
lookahead  distance  d  =  2  shiplengths. 
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Figure  7.   Vehicle  path  in  unstable  conditions  where  time 


constant  t, 


1.0  sec,  time  lag  T  =  1.0  sec,  and 


lookahead  distance  d  =  1.25  shiplength. 
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IV.    STABILITY  ENHANCEMENT 

A.   TWO -POINT  FORMULA 

The  time  lag  of  y  position  information  discussed  in 
chapter  II  has  thus  far  been  represented  by  a  single  piece  of 
information  used  in  the  control  law  with  some  delay  assigned 
to  it.  In  an  attempt  to  improve  stability  with  regard  to  this 
delay,  the  use  of  the  previous  two  positions  in  the  control 
was  examined. 

Equation  (2.15)  shows  the  representation  of  the  time 
delay,  and  equation  (2.16)  shows  its  effect  on  the  control 
law.  A  straight  line  approximation  is  applied  to  the  previous 
two  positions  as  shown  in  Figure  8.  The  general  y  position  is 
defined  by 

y(t)  =  axt+a2  (4.1) 

The   two  previous  positions  become 

y1=a1t1+a2  (4.2a) 

y2=axt2+a2  (4.2b) 

The   coefficients   ax   and  a2   are   stated  as 
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Figure  8.   Straight  line  approximation  using  two  previous 
positions . 
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at  =  il^  (4.3a) 


a.YA-y2t,  (4^3b) 


The  previous  time  values  are  represented  as 

t1  =  t-2T  (4.4a) 

t2  =  t-T  (4.4b) 

where  T  is  the  time  lag  associated  with  the  y  position 
information.  The  previous  positions,  in  terms  of  the  previous 
time  values  are 

y1=y(t-2T)  (4.5a) 

y7=y(C-T)  (4.5b) 


Substitution   and   algebra   gives   a   general   position 
expression  of 

y(t)  =2y(t-T)  -y{t-2T)  (4.6) 
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B.   STABILITY  ANALYSIS 

The  Laplace  transform  of  equation  (4.6)  is 

y(t)  =y(2e-Ts-e~2Ts)  (4.7) 

The  control  law  becomes 

b=k1\\i+k2v+k3r  +  ^y(2e-Ts-e-2Ts)  (4.8) 

and  the  linearized  equations  of  motion  become 

i|/=r  (4.9a) 


v  =  b1u2k1ty  +  (allu+b1u2k2)  v+  (a12u+b1u2k2)  r 
+bxu2^y(2e-Ts-e~2Ts) 


(4.9b) 


f  =  b2u2k1ty  +  (a12u+b2u2k2)  v+  (a22u+b2u2k2)  r 

k  (4    9c) 

+b2u2-^y(2e-Ts-e-2Ts) 


y=ui|f  +  v  (4.9d) 

These  motion  equations  reduce  to  the  state  space  form 

x  =  Ax  (4.10) 
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The  matrix  form  becomes 


v|  piU2^   axxu+bxu2k2   a12u^bxu2k3   bxu2-±  (2e-Ts-e~2Ts) 
b2u2k1   a21u+b2u2k2   a22u+b2u2k2    b2u2—±  {2e'Ts-e'2Ts) 


y\ 


The  characteristic  equation  of  the  form  [A-Is]=0  is 

s4+A3s3+^2s2^1s+(fl2s2+jb1s+fl0)D(2e"rs-e"2rs)  (4.11) 


where 


A3  =  -  (<311+<322^  u~  (b^+b^i)  u2 


A2  =  ia11a22-al2a21)  u2+  (a11b2-a21b1)  u3k3+  {a22bx-a12b2)  u3k2-b2u2k1 


Ax  =  -  (a21b1-a11b2)  u3kx 


Bx  =  -b2u2kx-  (a12b2-a22bx)  u3kx 
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B0  =   ^lA-^A)"^! 


H 


The  characteristic  equation  of  the  system  is  written  as 

l+KG(s)   =0  (4.12) 


where 


a(i).^^y|2.-"-e-»")  3 


s(s3+^>s2+^,s+^1) 


is  the  transfer  function,  and  we  have  denoted 

K=D  (4.14) 

With  s  =  ja>,  the  phase  angle  is  represented  by 

<t>  =  ZG(jco)  (4.15) 


where 

4>  =  Z(2e-^w-e-2r-7t0)  -Z(jo))  +l(-B2tt>2+B1jui+B0) 
-L  (-j(j33-A^2+A2j(j^+A1) 
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tan'1  -2sin(G)T')^sin(2a)T)  _n  +tan_1    B^   } 
2cos  (uT)  -cos  (2o)T)   2        Bn-B,G>2 

a  .  (4.16) 

.  .  -o)J-^2a)  . 

-tan"1  — 


The  Nyquist  stability  criterion  states  that  at  the  solution  of 

<j>(o>)  =-7i  (4.17) 

where  co  is  the  phase  cross  over  frequency  o>1#  the  gain  margin 
must  equal  1 

|JCG(jco1)  |=1  (4.18) 

where  K  =  D  =  1/d.   Equation  (4.17)  becomes 

.  -2sin(a)1D  +sin(2(o1D      n  .  ,      5,0),     v 

tan    -= 1 — —, 7^ — h~  " -?  ^tan"1  ( ±-±-  ) 

2cos  (w1D  -cos  (2co1D        2  S  -5  o) 


-tairM 


-wJ+^Wi 


■^q^+A, 


(4.19) 


Rearranged   as 


tan"1  ( ^_L_ )  -tan"1  ( 1      ^    1 


71  -tan"1  "2 

2  2cos(o)1D -cos(2<o1r) 


n  .  -2sin(a>1r) +sin(2o>,r) 

--tan"1- 
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Taking  the  tangent  of  both  sides  and  using  the  identity 


tan(x-y)=  tan U) -tan (y) 
1  +  tan (x)  tan(y) 


then  algebra  and  rearranging  gives 

B1<Ji1(A1-A3U)l)  -{BQ-B2u\)   (-cJi+A^) 

(B0-B2o>i)  {-A^I+AJ  +B1o1( -toJ+A^)  (4.20) 

2cos  (d)1r-cos  (2d)1D 

-2sin(co1r)  +sin(2a)1D 

The  stability  computations  require  an  initial  value  for 
co1.  This  is  accomplished  by  setting  the  time  lag  T=0,  and 
using 

(S0-S2Wi)  (-A3<i)?+A1)+B101(-<»J+A2ft>1)  =0       (4.21) 

rearranged  as 

(B^-Bi)  (tit*  (B1A2-B2A1-B0A2)  Wi+B0A1  =  0        (4.22) 

Setting  the  gain  margin  equal  to  1  gives 

1  K-BjwJ+B^+Bo)  (2e-ju'r-e-2ju'7)  |  _ 


J-jwJ-AjW^A-jo)^^ 


and  d  can  be  solved  for  by 
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(x>1\[(A1-A3u21)2+  (A2O1-0>i) 2 


In  order   to   facilitate   computation,    terms   are   grouped  as 
follows 


o2*fl1A2-S2A1-JB0A3 


i— "Ba 


Equation  (4.20)  becomes 


P1<o5  +  P2cj3+P3g)       2cos(cor)  -cos(2o)D 
a1o)4+a2a)2+a3        -2sin  (a>r)  +  sin  (2<i>D 


Put    in   the    form 
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(4.24) 


f(to)  =0  (4.25) 


equation  (4.24)  becomes 

[  (p1co5  +  P2co3  +  P3io)  (-2sin(o>T)  +sin(2o>r)  ) 


(4.26) 
[  (a1o)4+a2a)z+a3)  (2cos  (uT)  -cos  (2uD  )  ]  =0 


Newton's  iteration  is  used  as 


f(co,  .) 

«*  =  ^-i-  17;     \  (4-27) 


Where    the    function   of    co   is 

f(co)  =  (p1(o5  +  p2o)3  +  p3o))  (-2sin(oD  ) 


(P1co5  +  P2co3  +  P3g>)  (sin(2wD  )  (4    28) 

(a1a)4+a2G)2+a3)  (2cos(oD) 
(a1a)4+a2a)2+a3)  (cos(2coD) 


and  the  derivative  function  is 


f"'(co)  =  (5P1o)4+3P2&)2  +  P3)  (-2sin(o>r)  ) 
-(p1o)5  +  P2(i)3  +  P1a))  (2Tcos(uT) 
+  (5p1o)4+3p2a)2  +  P1)  (sin(2«D) 
+  (P1co5  +  P2to3+P3io)  (2Tcos(2v>T)  )  (4    2£) 

+  (4a1o3+2a2G))  (-2cos  ((i>D  ) 
+  (a1to4+a2a)2+a3)  (2rsin(coD) 
+  (4a1to3+2a2co)  (cos(2o)D) 
-  (a1o)4+a2o)2+a3)  (2Tsin  (2 cor)  ) 
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Appendix  C  presents  the  Fortran  program  used  to  determine 
stability  characteristics  as  a  function  of  system  time 
constant  and  time  lag. 

C.   RESULTS  AND  DISCUSSION 

Values  for  the  minimum  lookahead  distance  required  for 
controller  stability  have  been  computed  as  a  function  of  the 
system  time  constant  and  the  y  position  information  time  lag 
using  the  two-point  formula.  Figures  9,  10,  and  11  show 
critical  lookahead  distance  plotted  against  time  constant  at 
three  values  of  time  lag  for  both  the  two-point  method  as  well 
as  the  single  point  (Nyquist)  method  discussed  in  chapter  II. 
These  figures  show  an  overall  decrease  in  required  lookahead 
distance  as  a  result  of  the  use  of  the  two-point  memory  model. 
An  exception  to  this  occurs  at  low  values  of  tc.  However, 
these  values  are  to  be  avoided  in  practice,  as  explained  in 
the  previous  chapter. 

Figure  12  through  15  show  the  critical  lookahead  distance 
plotted  against  time  lag  for  different  values  of  the  time 
constant  using  both  the  two-point  and  single  point  methods. 
Again  the  two-point  method  produces  superior  results  to  the 
single  point  method,  with  the  exception  of  small  time 
constants  and  large  time  lags.  The  same  information  is 
summarized  in  Figure  16  through  19,  where  we  show  the  ratio 
d2/d1  (critical  lookahead  distance  using  the  two-point  method 
divided  by  the  critical  lookahead  distance  using  the  single 
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point  method)  verses  time  constant  and  time  lag.  It  can  be 
seen  that  the  use  of  the  two-point  method  can  result  in 
improvement  in  the  region  of  stability  of  over  20%. 
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2.2 


Figure  9.   Critical  lookahead  distance  d  versus  time  constant 


with   time   lag   T 


0.5   sec   (two-point 


method=solid,  single  point=dashed) . 
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0.6 


Figure  10.  Critical  lookahead  distance  d  versus  time  constant 
tc  with  time  lag  T  =  1.0  sec  (two -point 
method=solid#  single  point=dashed) . 
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0.8 


Figure  11.  Critical  lookahead  distance  d  versus  time  constant 
tc  with  time  lag  T  =  1.5  sec  (two-point 
method=solid,  single  point=dashed) . 
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Figure  12.  Critical  lookahead  distance  d  versus  time  lag  T 
with  time  constant  tc  =  0.8  sec  (two-point 
method=solid,  single  point=  dashed) . 
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2.1) 


Figure  13.  Critical  lookahead  distance  d  versus  time  lag  T 
with  time  constant  tc  =  1.0  sec  (two-point 
method=solid#  single  point=  dashed) . 
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Figure  14.  Critical  lookahead  distance  d  versus  time  lag  T 
with  time  constant  tc  =  1.5  sec  (two -point 
method=solid,  single  point=dashed) . 
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Figure  15.  Critical  lookahead  distance  d  versus  time  lag  T 
with  time  constant  tc  =  2.5  sec  (two-point 
method=solid,  single  point=dashed) . 
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1.2 


Figure  16.  Ratio  d^/d^  versus  time  constant  tc  with  time 
lag  =  1.0  sec. 
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Figure  17.  Ratio  d2/d1  versus  time  constant  tc  with  time 
lag  =  1.5  sec. 
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figure  18.  Ratio  d2/d1  versus  time  lag  T  with  time  constant  = 
1.0  sec 
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Figure  19.  Ratio   d2/d1  versus  time  lag  T  with  time 
constant  =  1.5  sec. 
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V.    CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

Nonlinear  bifurcation  theory  has  been  applied  to  the  NPS 
autonomous  underwater  vehicle  through  a  Hopf  bifurcation 
analysis  of  the  vehicle's  navigation  guidance  and  control 
scheme.  It  has  been  shown  that  under  normal  conditions  the 
vehicle  will  operate  in  the  supercritical  region,  and 
therefore  the  unpredictable  behavior  associated  with  the 
subcritical  case  should  not  to  be  expected. 

Secondly,  the  control  scheme  was  modified  by  the  use  of  a 
single  stage  memory  model  which  incorporates  the  two  previous 
vehicle  positions  in  the  linearized  control  law.  It  has  been 
shown  that  the  use  of  this  memory  model  reduces  the  vehicle's 
required  lookhead  distance  for  control  stability  within  the 
normal  operating  range. 

B .  RECOMMENDATIONS 

In  the  event  that  not  all  of  the  state  information 
required  is  directly  available  through  measurements,  the 
control/guidance  scheme  would  include  an  observer.  The  effect 
of  the  observer  dynamics  on  the  higher  order  nonlinear  terms, 
and  therefore  on  the  bifurcation  limit  cycle  stability,  should 
be  examined. 
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APPENDIX  A 


PROGRAM  HOPF.FOR 

Hopf  Bifurcation  Analysis 

Critical  Point:  Exact 

Third  Order  Expansions:   First  Order  Approximation 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , K4 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
&  NRDOT,NVDOT,KlP,K2P 

DOUBLE  PRECISION  Ml 1 , M12 , Ml  3 , Ml  4 , M21 , M22 , M23 , M24 , 

1  M31,M32,M33,M34,M41,M42,M4  3,M4  4, 

2  N11,N12,N13,N14,N21,N2  2,N2  3,N2  4, 

3  N31,N32,N33,N34,N41,N4  2,N4  3,N4  4, 

4  L21fL2  2,L2  3,L2  4,L31,L32,L33,L34, 

5  L41,L42,L43,L44 

DIMENSION  A(4,4) ,T(4,4) ,TINV(4,4) , FVl ( 4 ) , I VI ( 4 ) , ZZZ ( 4 , 4 ) 
DIMENSION  WR( 4 ) ,WI ( 4 ) ,TSAVE( 4,4) , TLUD( 4,4  )  , IVLUD( 4 ) ,SVLUD( 4 ) 
DIMENSION  ASAVE(4,4) 

OPEN  ( ll,FILE='HOPF.DAT' , STATUS- 'NEW ) 
OPEN  ( 12,FILE='PERI .DAT' , STATUS- ' NEW ) 
OPEN  ( 1 3, F I LE=' ALPH10.DAT' , STATUS- 'NEW ) 

Vehicle  Parameters 

WEIGHT-435.0 
IZ     =45.0 
L      =7.3 
RHO    =1.94 
G      =32.2 
XG     =0.0104 
MASS   -WEIGHT/G 
U      =5.0 
RATIO  =0.0 
DO     =0.4 

WRITE  (*,1006) 

READ   ( *  ,  *  )  DO 

YRDOT  =-0. 00178*0. 5*RHO*L**4 
YVDOT  —  0.03430*0. 5*RHO*L**3 
YR  =+0.01187*0. 5*RHO*L**3 
YV  =-0.03896*0. 5*RHO*L**2 
YDRS  =+0. 02345*0. 5*RHO*L**2 
YDRB  =+0. 02345*0. 5*RHO*L**2 
NRDOT  =-0. 00047*0. 5*RHO*L**5 
NVDOT  =-0. 00178*0. 5*RHO*L**4 
NR  =-0. 01022*0. 5*RHO*L**4 
NV     =-0. 00769*0. 5*RHO*L**3 
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*)R5        —  U.JJ/*U.UZJ4:>*U.  b*RHU*L'"fJ 

*>RB      -+0.283*0. 02345*0. 5*RH0*L** 3 

.1      «(  IZ-NRDOT)* (MASS-YVDOT) - 

( MASS*XG-YRDOT ) * ( MASS*XG-NVDOT ) 
kll-( ( IZ-NRDOT )*YV-(MASS*XG-YRDOT)*NV)/DH 
i\12=(  (  IZ-NRDOT) *(-MASS+YR)- 

(MASS*XG-YRDOT)*(-MASS*XG+NR) )/DH 
!v21=( (MASS-YVDOT) *NV-( MASS *XG-NVDOT) *YV)/DH 
[^22=  (  (  MASS-YVDOT  )  *  (  -MASS*XG+NR  )  - 

(MASS*XG-NVDOT)*(-MASS+YR) )/DH 
•Jll=(  (  IZ-NRDOT)  * YDRS- ( MASS*XG-YRDOT )  *NDRS  )/DH 
iil2=(  (  IZ-NRDOT  )*YDRB-(MASS*XG-YRDOT)*NDRB)/DH 
■J21-  (  (  MASS-YVDOT  )  *NDRS-  (  MASS*XG-NVDOT  )  *  YDRS  )  /DH 
:l22-(  (MASS-YVDOT) *NDRB-(MASS*XG-NVDOT) *YDRB )/DH 


Jl»BBll+RATIO*BBl2 
J2=BB21+RATI0*BB22 


UTE  (*,1005) 

SAD  (  *  ,  *  )  TL 

UTE  (*,1001) 

!AD  (  *  ,  *  )  TC 


■-1.D-10 

l?MAX=10  00 

[•L=0.5 

)0    50    M-1,100 

!JD=TL*L/U 

?C=0.5 

:)    40    K-1,100 

TCD=TC*L/U 

POLEC=l .0/TCD 

AD1=3.0*POLEC 

AD2=3.0*POLEC**2 

AD3=POLEC**3 

A1=BB1*U*U 

B1=BB2*U*U 

Cl— ADl-(AAll+AA22)*U 

A2=(BB1*AA2  2-BB2*AA12)*U**3 

B2=(BB2*AAll-BBl*AA21 )*U**3 

Kl=AD3/( (BB2*AA11-BB1*AA21)*U**3) 

C2=AD2-(AA11*AA2  2-AA12*AA21)*U**2+BB2*U*U*K1 

K2=(C1*B2-C2*B1)/(A1*B2-A2*B1) 

K3=(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 

Al=  (AA11*BB2-AA21*BB1 )*U*U*U*K1 

A2=  (AA11*AA2  2-AA12*AA21 ) *U*U+ ( AAll *BB2-AA21 *BB1 )*U*U*U*K3- 

(AA22*BB1-AA12*BB2)*U*U*U*K2-BB2*U*U*K1 
A3— (AAll+AA22)*U-(BBl*K2+BB2*K3)*U*U 
B0=  (AA11*BB2-AA21*BB1)*U*U*U*U*K1 
B1=-(AA12*BB2-AA22*BB1)*U*U*U*K1-BB2*U*U*U*K1 
B2=-BB1*U*U*K1 

Compute  Initial  Approximation  For  Omega  (TL=0) 
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AQ=B2*A3-B1 

BQ=Bl*A2-B2*Al-B0*A3 

CQ=BO*Al 

DET=BQ**2-4 .DO*AQ*CQ 

IF  (DET.LT.O.DO)  GO  TO  40 

ZT1=(-BQ+DSQRT(DET) )/(2.0D0*AQ) 

ZT2=(-BQ-DSQRT(DET) )/(2.0D0*AQ) 

IF  (ZT1.GE.0.D0)  OM=DSQRT(ZTl) 

IF  (ZT2.GE.0.D0)  OM=DSQRT( ZT2 ) 

ALPHA1-AQ 

ALPHA2=BQ 

ALPHA3=CQ 

BETAl— B2 

BETA2*B0+B2*A2-Bl*A3 

BETA3=Bl*Al-B0*A2 
r 
C       Compute  Exact  Omega  Using  Newton's  Iteration 

9  3    F=(BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD)*DSIN(OMOLD*TLD) 
&    +(ALPHAl*OMOLD**4+ALPHA2*OMOLD**2+ALPHA3)*DCOS(OMOLD*TLD) 

FPRIME=(5  D0*BETAl*OMOLD**4+3.D0*BETA2*OMOLD**2+BETA3) 
&    *DSIN(OMOLD*TLD)+(BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD) 
&    *TLD*DCOS(OMOLD*TLD)+(4.D0*ALPHAl*OMOLD**3+2.D0*ALPHA2*OMOLD 
&    *DCOS(OMOLD*TLD)-(ALPHAl*OMOLD**4+ALPHA2*OMOLD**2+ALPHA3) 

£,    *TLD*DSIN(OMOLD*TLD) 

IF  (FPRIME.EQ.0.D0)  STOP  1112 
IF  (F.EQ.0.D0)  GO  TO  92 
OMNEW=OMOLD-F/FPRIME 
OMDI F=DABS ( OMNEW-OMOLD ) 
IF  (OMDIF.LE.EPS)  GO  TO  92 
OMOLD=OMNEW 
IT=IT+1 

IF  ( IT.GT. ITMAX)  STOP  1111 
GO  TO  9  3 
C 

92    OM=OMNEW 

XDNUM=DSQRT(  ( B0-B2 *OM* *  2 ) *  *  2+ ( Bl *OM ) *  *  2 ) 
XDDEN=OM*DSQRT(  ( Al-A3*OM* *  2 ) *  *  2+ ( A2 *OM-OM* *  3 ) *  *  2 ) 
XD=XDNUM/XDDEN 

C  ,    • 

C       Start  Hopf  Bifurcation  Analysis 

C 

C       Evaluate  Nonlinear  Rudder  Expansion  Coefficients 

C 


C 


KlP=Kl-Kl*TLD*U/XD 
K2P=K2-K1*TLD/XD 

DPPV=-  (l.D0/(3.D0*D0**2) ) *3 . D0*KlP*KlP*K2P 

&        +  0.5D0*K1*TLD/XD  +  3 . D0*Kl * ( TLD*U ) *  * 3/( 3 . D0*XD* *  3 
DPVV=-  (l.D0/(3.D0*D0**2)  ) *  3 . D0*KlP*K2P*K2P 

&        +  3.D0*U*TLD*TLD*TLD*Kl/( 3.D0*XD**3) 
DPPR=-  (l.D0/(3.D0*D0**2) ) *3 . D0*KlP*KlP*K3 
DPRR=-  (l.D0/(3.D0*D0**2) ) *  3 . D0*KlP*K3*K3 
DPPY=-  (l.D0/(3.D0*D0**2)  ) *  3 . D0*KlP*KlP*Kl/XD 

&        -  3.D0*TLD*TLD*U*U*Kl/(3.D0*XD**3) 

DPYY=-  (l.D0/(3.D0*D0**2)  ) *  3 . D0*KlP*Kl *Kl/( XD* *2 ) 

&        +  3.D0*TLD*U*Kl/(3.D0*XD**3) 

DVVR=-  (l.D0/(3.D0*D0**2)  ) *  3 . D0*K2P*K2P*K3 
DVRR=-  (l.D0/(3.D0*D0**2)  ) *  3 . D0*K2P*K3*K3 
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-  3.D0*K1*TLD*TLD/(3.D0*XD**3) 

DVYY=-  (l.D0/(3.D0*D0**2) ) *3 . D0*Kl*Kl*K2P/( XD* *2 ) 

+  3.D0*TLD*Kl/( 3.D0*XD**3) 
DRRY=-  (l.DO/( 3.D0*D0**2)  ) *  3 . DO*Kl *K3*K3/XD 
DRYY—  (l.D0/(3.D0*D0**2)  ) *  3 . DO*Kl*Kl *K3/( XD**2 ) 
DPVR=-  (l.DO/( 3.D0*D0**2) ) *6 . D0*KlP*K2P*K3 
DPVY=-  ( l.DO/( 3.D0*D0**2) ) *6 . D0*KlP*Kl*K2P/XD 

-  6.D0*TLD*TLD*U*Kl/( 3.D0*XD**3 ) 
DPRY=-  (1 .D0/( 3.D0*D0**2) ) *6 . DO*KlP*Kl *K3/XD 
DVRY=-  (l.DO/( 3.D0*D0**2) )  *6 . DO*Kl *K2P*K3/XD 
DPPP=-  (l.DO/( 3.D0*D0**2))*l.D0*KlP*KlP*KlP 

+  Kl*TLD*U/(6.D0*XD)  +  ( Kl * ( TLD*U ) *  *  3 ) /( 3 . DO *XD* *  3 
DVW=-  (l.D0/(3.D0*D0**2)  )  *1  .  D0*K2P*K2P*K2P 

+  Kl*(TLD**3)/( 3.D0*XD**3) 
DRRR=-  (l.D0/(3.D0*D0**2) )*1.D0*K3*K3*K3 
DYYY=-  (l.D0/(3.D0*D0**2)  ) * ( Kl/XD ) *  * 3-Kl/( 3.D0*XD**3) 

-  Kl/(3.D0*XD**3) 

Evaluate  Transformation  Matrix  of  Eigenvectors 

PSI  =  0.D0 

Y  =0.D0 

V  =0.D0 
DPSI=K1P 
DV   =  K2P 
DR   =K3 

DY   =K1*XD/(XD**2+Y**2) 

|(lf  D-0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2,l)=        BB1*U*U*DPSI 

A(2,2)=AA11*U+BB1*U*U*DV 

A(2, 3)=AA12*U+BB1*U*U*DR 

A(2,4)=        BB1*U*U*DY 

A(3fl)=        BB2*U*U*DPSI 

A(  3,2)=AA21*U+BB2*U*U*DV 

A(  3, 3)=AA22*U+BB2*U*U*DR 

A(3f4)=        BB2*U*U*DY 

A(4, 1 )=U*DCOS(PSI )-V*DSIN(PSI ) 

A(4,2)=   DCOS(PSI) 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

DO  11  1=1,4 

DO  12  J=l,4 

ASAVE( I, J)=A( I, J) 

CONTINUE 
CONTINUE 

CALL  RG(4,4,A,WR,WI , 1 , ZZZ , IVl , FVl , I  ERR ) 
CALL  DSOMEG( IEV,WR,WI , OMEGA , CHECK ) 
OMEGA0=OMEGA 
DO  5  1=1,4 

T(I,1)=  ZZZ( I,IEV) 

T(I ,2)=-ZZZ( I , IEV+1 ) 
CONTINUE 

IF  (IEV.EQ.l)  GO  TO  13 
IF  (IEV.EQ.2)  GO  TO  14 
IF  (IEV.EQ.3)  GO  TO  15 
STOP  3004 
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15 


13 


16 
17 


T(I,3)=ZZZ( 1,1) 

T(I,4)=ZZZ(I,4) 
CONTINUE 
GO  TO  17 
DO  7  1=1,4 

T(I,3)=ZZZ(I,1) 

T(I  ,4)=ZZZ(I,2) 
CONTINUE 
GO  TO  17 
DO  16  1=1,4 

T( I,3)=ZZZ(I,3) 

T( 1,4 )=ZZZ( 1,4) 
CONTINUE 
CONTINUE 

Normalization  of  the  Critical  Eigenvector 

CALL  NORMAL(T) 

Definition  of  Mij 


M11=T(1, 

,1) 

M21=T(2, 

.1) 

M31=T(3, 

.1) 

M41=T(4, 

,1) 

M12=T(1, 

,2) 

M22=T(2, 

-2) 

M32=T(3( 

-2) 

M42=T( 4, 

-2) 

M13=T(1, 

-3) 

M23=T(2, 

.3) 

M33=T(3, 

-3) 

M43=T(4( 

,3) 

M14=T(1, 

,4) 

M24=T(2( 

,4) 

M34=T(3, 

,4) 

M44=T(4( 

r4) 

Definition  of  Lij 

L21=  DPPV*M11*M11*M21 

&  +  DPRR*M11*M31*M31 

&  +  DVVR*M21*M21*M31 

&  +  DVYY*M21*M41*M41 

&  +  DPVR*Mll*M21*M31 

&  +  DVRY*M21*M41*M31 

&  +  DRRR*M31*M31*M31 


DPVV*M11*M21*M21 
DPPY*M11*M11*M41 
DVRR*M21*M31*M31 
DRRY*M31*M31*M41 
DPVY*M11*M21*M41 
DPPP*Mll*Mll*Mll 
DYYY*M41*M41*M41 


DPPR*Mll*Mll*M31 
DPYY*M11*M41*M41 
DVVY*M21*M21*M41 
DRYY*M31*M41*M41 
DPRY*M11*M41*M31 
DVVV*M21*M21*M21 


PPV 
PVV 
PPR 
PRR 
PPY 
PYY 
VVR 
VRR 
VVY 
VYY 
RRY 


M11*M11 
M12*M21 
Mll*Mll 
M12*M31 
Mll*Mll 
M41*M41 
M21*M21 
M22*M31 
M21*M21 
M22*M41 
M31*M31 


*M22 
*M21 
*M32 
*M31 
*M42 
*M12 
*M32 
*M31 
*M42 
*M41 
*M42 


2.0*M11 
2.0*M11 
2.0*M11 
2.0*Mll 
2.0*M11 
0*Mll 
0*M31 
0*M31 
0*M41 
0*M41 
0*M41 


*M12*M21 
*M21*M22 
*M12*M31 
*M31*M32 
*M12*M41 
*M41*M42 
*M21*M22 
*N32*M21 
*M21*M22 
*M42*M21 
*M31*M32 
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RYY  -  M32*M41*M41  +  2 . 0*M4 1 *M42*M31 
PVR  -  M11*M21*M32+M31*(M11*M22+M12*M21) 
PVY  -  M11*M21*M42+M41*(M11*M22+M12*M21) 
PRY  -  M11*M41*M32+M31*(M11*M42+M12*M41 ) 
VRY  =  M21*M41*M32+M31*(M21*M42+M22*M41 ) 
PPP  -  3.0*Mll*Mll*Ml2 
VVV  «  3.0*M21*M21*M22 
RRR  -  3.0*M31*M31*M32 
YYY  -  3.0*M41*M41*M42 

L22=DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 
+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 
+DRRR*RRR+DYYY*YYY 


PPV 
PVV 
PPR 
PRR 
PPY 
PYY 
WR 
VRR 
VVY 
VYY 
RRY 
RYY 
PVR 
PVY 
PRY 
VRY 
PPP 
VVV 
RRR 
YYY 


M12*M1 
Mll*M2 
M12*M1 
Mll*M3 
M12*M1 
Mll*M4 
M22*M2 
M21*M3 
M22*M2 
M21*M4 
M32*M3 
M31*M4 
M12*M2 
M12*M2 
M12*M4 
M22*M4 
3.0*Ml 
3.0*M2 
3.0*M3 
3.0*M4 


2*M21 
2*M22 
2*M31 
2*M32 
2*M41 
2*M42 
2*M31 
2*M32 
2*M41 
2*M42 
2*M41 
2*M42 
2*M31 
2*M41 
2*M31 
2*M31 
1*M12 
1*M22 
1*M32 
1*M42 


2.0*M11 
2.0*M12 


0*M11 

0*M12 
0*M11 
0*M12 
0*M21 
0*M22 
0*M21 
0*M22 
+    2.0*M31 
+    2.0*M32 
+    M32*(M1 
+    M42*(M1 
+    M32*(M1 
+    M32*(M2 
*M12 
*M22 
*M32 
*M42 


*M12*M22 

*M21*M22 

*M12*M32 

*M31*M32 

*M12*M42 

*M41*M42 

*M22*M32 

*M31*M32 

*M22*M42 

*M41*M42 

*M32*M42 

*M41*M42 

1*M22+M12*M21 ) 

1*M22+M12*M21 ) 

1*M42+M12*M41 ) 

1*M42+M22*M41 ) 


L2  3=DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 
+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 
+DRRR*RRR+DYYY*YYY 


L24=  DPPV*Ml2*Ml2*M22 

+  DPRR*M12*M32*M32 

+  DVVR*M22*M22*M32 

+  DVYY*M22*M42*M42 

+  DPVR*M12*M22*M32 

+  DVRY*M22*M42*M32 


DPVV*M12*M22*M22 
DPPY*M12*M12*M42 
DVRR*M22*M32*M32 
DRRY*M32*M32*M42 
DPVY*M12*M22*M42 
DPPP*M12*M12*M12 


+  DRRR*M32*M32*M32  +  DYYY*M42*M42*M42 

L31=L21 
L32=L22 
L33=L23 
L34=L24 

L21«L21*BB1*U*U 
L22=L22*BBl*U*U 
L23=L23*BBl*U*U 
L24=L24*BBl*U*U 
L31=L31*BB2*U*U 
L32=L32*BB2*U*U 


DPPR*M12*M12*M32 
DPYY*M12*M42*M42 
DVVY*M22*M22*M42 
DRYY*M32*M42*M42 
DPRY*M12*M42*M32 
DVVV*M22*M22*M22 
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b33  =  LJ0"DDZ''U"U 

L34=L34*BB2*U*U 

L41=-0.5*Mll*Mll*(M21+U*Mll/3.0) 
L42=-Mll*(Ml2*M21+0.5*Mll*M22+0.5*U*Ml2*Ml2) 
L4  3=-M12*(M11*M2  2+0.5*M12*M21+0.5*U*M11*M12) 
L4  4— 0.5*Ml2*Ml2*(M22+U*Ml2/3.0) 

Invert  Transformation  Matrix 

DO  2  1=1,4 
DO  3  J=l,4 

TINV(I,J)«0.0 

TSAVE( I , J)=T( I , J) 

3  CONTINUE 
2    CONTINUE 

CALL  DLUD (4,4, TSAVE , 4 , TLUD , I VLUD ) 
DO  4  1=1,4 

IF  (IVLUD(I) .EQ.O)  STOP  3003 

4  CONTINUE 

CALL  DILU( 4, 4, TLUD, IVLUD,SVLUD) 
DO  8  1=1,4 
DO  9  J=l,4 

TINV( I , J)=TLUD( I, J) 
9      CONTINUE 
8    CONTINUE 

Check  Inv(T)*A*T 

IMULT=2 

IF  (IMULT.EQ.l)  CALL  MULT ( TINV , ASAVE , T ) 

IF  (IMULT.EQ.O)  STOP 

Definition  of  Nij 

Nll=TINV(l,l ) 
N21=TINV(2,1) 
N31=TINV( 3,1 ) 
N41=TINV( 4,1 ) 
N12=TINV( 1,2) 
N22=TINV( 2,2 ) 
N32=TINV( 3,2) 
N42=TINV( 4,2) 
N13=TINV(1,3) 
N23=TINV(2, 3) 
N33=TINV( 3,3) 
N43=TINV( 4,3) 
N14=TINV( 1,4 ) 
N24=TINV(2,4 ) 
N34=TINV( 3,4) 
N44=TINV( 4,4 ) 

R11=N12*L21+N13*L31+N14*L41 
R12=N12*L2  2+N13*L32+N14*L4  2 
R13=N12*L2  3+N13*L3  3+N14*L4  3 
R14=N12*L2  4+N13*L3  4+N14*L4  4 
R21=N22*L21+N2  3*L31+N2  4*L41 
R22=N22*L22+N2  3*L32+N2  4*L42 
R2  3=N2  2*L2  3+N2  3*L3  3+N2  4*L4  3 
R24=N22*L2  4+N2  3*L3  4+N2  4*L4  4 
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Evaluate  Aipna-  ana 

umega • 

DINC=0.001 

XDR 

-XD+DINC 

XDL 

-XD-DINC 

XD 

«XDR 

PSI 

=  0.D0 

Y 

-0.D0 

V 

-O.DO 

DPSI=K1P 

DV 

=  K2P 

DR 

=  K3 

DY 

=Kl*XD/(XD**2+Y**2) 

A(l 

D-0.0D0 

A(l 

2)=0.0D0 

A(l 

3)-1.0D0 

A(l 

4)=0.0D0 

A(2 

1)=         BB1*U*U*DPSI 

A(2 

2)=AAll*U+BBl*U*U*DV 

A(2 

3)=AAl2*U+BBl*U*U*DR 

A(2,4)  =        BB1*U*U*DY 

A(3,l)  =        BB2*U*U*DPSI 

A(  3,2)=AA21*U+BB2*U*U*DV 

A( 3, 3)=AA22*U+BB2*U*U*DR 

A(3,4)  =        BB2*U*U*DY 

A(4f 1 )=U*DCOS(PSI )-V*DSIN(PSI ) 

A(4,2)  =   DCOS(PSI) 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

CALL  RG( 4,4,A,WR,WI , 0,ZZZ, IVl , FVl , I  ERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
ALPHR=DEOS 
OMEGR=FREQ 

XD=XDL 
PSI  =0.D0 

Y  =0.D0 

V  =0.D0 
DPSI=K1P 
DV   =K2P 
DR   =K3 

DY   =K1*XD/(XD**2+Y**2) 

A(1,1)=0.0D0 

A(1,2)-0.0D0 

A(1,3)=1.0D0 

A(1,4)-0.0D0 

A(2,l)=        BB1*U*U*DPSI 

A(2,2)=AAll*U+BBl*U*U*DV 

A(2, 3)=AAl2*U+BBl*U*U*DR 

A(2,4)=        BBl*U*U*DY 

M3,l)  =        BB2*U*U*DPSI 

M3,2)=AA21*U+BB2*U*U*DV 

M  3 , 3 )=AA22*U+BB2*U*U*DR 

M3,4)=        BB2*U*U*DY 

M4,l)=U*DCOS(PSI  )-V*DSIN(PSI  ) 

M4,2)  =   DCOS(PSI) 

■  4,3)«0.0D0 

14,  4)=0.0D0 
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CALL  DSTABL(DEOS,WR,WI , FREQ ) 

ALPHL=DEOS 

OMEGL-FREQ 

DALPHA=(ALPHR-ALPHL)/(XDR-XDL) 
DOMEGA= ( OMEGR-OMEGL ) / ( XDR-XDL ) 

Evaluation  of  Hopf  Bifurcation  Coefficients 

COEF1=3.0*R11+R13+R22+3.0*R2  4 

COEF2=3.0*R21+R2  3-R12-3.0*R14 

AMU2  =-COEFl/(8.0*DALPHA) 

BETA2=0.25*COEF1 

TAU2  =- ( COEF2-DOMEGA*COEFl/DALPHA)/( 8 . 0*OMEGA0 ) 

PERD   =2.0*3.1415927/OMEGA0 

PER   =PERD*U/L 

X=XD/L 

K4=Kl/XD 

WRITE  (11,*)  TC,COEFl 


1 

WRITE  (12, 

*)  TC,PER 

1 

WRITE  (13, 

*)  TC,DALPHA 

TC=TC+0.015 

40  CONTINUE 

c 

TL=TL+0. 02 

£ 

50  CONTINUE 

'  1001  FORMAT  ( '  ENTER  TIME  CONSTANT' ) 

1005  FORMAT  ('  ENTER  TIME  LAG') 

1006  FORMAT  ('  ENTER  DO') 
-  2002  FORMAT  (5E15.5) 

:  END 

SUBROUTINE  DSOMEG( IJK,WR,WI , OMEGA , CHECK ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 
CHECK=-1.0D+25 
DO  1  1=1,4 

IF  (WR( I ) .LT. CHECK)  GO  TO  1 

CHECK=WR( I ) 

IJ  =  I 
1  CONTINUE 

OMEGA=DABS(WI(IJ) ) 

IF  (WI(IJ) .GT.0.D0)  IJK=IJ 

IF  (WI(IJ) .LT.0.D0)  IJK=IJ-1 

RETURN 

END 


SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4) ,WI(4 ) 
DEOS=-l .OD+20 
DO  1  1=1,4 

IF  (WR(I ) .LT.DEOS)  GO  TO  1 

DEOS=WR(I ) 

IJ  =  I 
CONTINUE 
OMEGA=WI( I  J) 
OMEGA=DABS ( OMEGA ) 
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ND 

UBROUTINE  NORMAL (T) 

MPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

IMENSION  T( 4,4) ,TNOR(4,4) 

FAC=T(1,1 )**2+T(l,2)**2 

F  (DABS(CFAC) .LE. (l.D-10) )  STOP  4001 

NOR(l,l)=l.D0 

NOR (2,1)-(T(1,1)*T(2,1)+T(2,2)*T(1,2)) /CFAC 

'NOR (3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2)) /C  FAC 

NOR(4,l)=(T(l,l)*T(4,l)+T(4,2)*T(l,2) )/CFAC 

WOR(l,2)-0.D0 

'NOR(2,2)=(T(2,2)*T(l, 1 )-T( 2 , 1 ) *T( 1 , 2 ) )/CFAC 

'NOR  ( 3 , 2 ) » ( T ( 3 , 2 ) *  T ( 1 , 1 ) -T ( 3 , 1 ) *T ( 1 , 2 ) ) /CFAC 

'NOR (4,2)=(T(4,2)*T(1,1)-T(4,1)*T(1,2))/C  FAC 

>0  1  1=1,4 

DO  2  J-1,2 

T( I , J)=TNOR( I , J) 

CONTINUE 
ONTINUE 
;ETURN 
IND 


UBROUTINE  MULT(TINV,A,T) 
MPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
•IMENSION  TINV(4,4),A(4,4),T(4,4),Al(4,4),A2(4,4 
»0  1  1  =  1,4 
DO  2  J-1,4 

A1(I, J)=0.D0 
A2( I, J)=0.D0 
CONTINUE 
'ONTINUE 
»0  3  1  =  1,4 
DO  4  J-1,4 
DO  5  K=l,4 

Al( I , J)=A( I , K)*T(K, J)+Al( I , J) 
CONTINUE 
CONTINUE 
lONTINUE 
>0  6  1  =  1,4 
DO  7  J-1,4 
DO  8  K-1,4 

A2( I , J)=TINV( I ,K)*A1(K, J)+A2( I , J) 
CONTINUE 
CONTINUE 
JONTINUE 
)0  11  1  =  1,4 

WRITE  (*,101)  (A( I , J) , J-1,4) 
:ONTINUE 
)0  12  1  =  1,4 

WRITE  (*,101)  (T( I, J) , J=l,4 ) 
:ONTINUE 
>0  10  1  =  1,4 

WRITE  (*,101)  (A2(I,J) , J-1,4) 
ONTINUE 

WRITE  (*,101)  A2(l,l) 
IETURN 

?ORMAT  (4E15.5) 
3ND 
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APPENDIX  B 


PROGRAM  SIMU.FOR  (VEHICLE  PATH  SIMULATION) 
REAL  Kl , K2, K3 , L , NR , NV , NDRS , NDRB , I Z, MASS, 
&      NRDOT,NVDOT,KlP,K2P 

OPEN  (2  0,FILE='SIMl .DAT' , STATUS- 'NEW ) 
OPEN  ( 3 0,FILE=' SIM2.DAT' , STATUS- 'NEW ) 

WEIGHT=435.0 

IZ     =45.0 

L      =7.3 

RHO    =1.94 

G      =32.2 

XG     =0.0104 

MASS   =WEIGHT/G 

U      =5.0 

RATIO  =0.0 

DO     =0.4 

XG     =XG/L 

MASS   =MASS/( 0. 5*RHO*L**3) 

IZ     =IZ/( 0. 5*RHO*L**5) 

YRDOT  =-0.00178 

YVDOT  =-0.0  34  3  0 

YR     =+0.01187 

YV     =-0.03896 

YDRS   =+0.02345 

YDRB   =+0.02345 

NRDOT  =-0.00047 

NVDOT  =-0.00178 

NR     =-0.01022 

NV     =-0.00769 

NDRS   =-0.337*0.02345 

NDRB   =+0.283*0.02345 

DH   =  (  IZ-NRDOT)*(MASS-YVDOT)- 

&       ( MASS*XG-YRDOT ) * ( MASS*XG-NVDOT ) 
AA11=( ( IZ-NRDOT)*YV-(MASS*XG-YRDOT)*NV)/DH 
AA12=( ( IZ-NRDOT) *(-MASS+YR )- 

&       (MASS*XG-YRDOT)*(-MASS*XG+NR) )/DH 
AA21= ( ( MASS-YVDOT ) *NV- ( MASS*XG-NVDOT ) * YV ) /DH 
AA22= ( ( MASS-YVDOT ) * ( -MASS*XG+NR ) - 

&       (MASS*XG-NVDOT)*(-MASS+YR) )/DH 
BB11=( ( IZ-NRDOT) * YDRS- ( MASS* XG- YRDOT )* NDRS ) /DH 
BB12=( ( IZ-NRDOT) *YDRB-(MASS*XG-YRDOT) *NDRB )/DH 
BB21=( (MASS- YVDOT )*NDRS-( MASS *XG-NVDOT)* YDRS )/DH 
BB22= ( ( MASS-YVDOT ) *NDRB- ( MASS*XG-NVDOT ) * YDRB )/DH 

BB1=BB11+RATI0*BB12 
BB2=BB21+RATIO*BB2  2 

WRITE  (*,*)  'ENTER  TC,TL,D' 

66 


ILE-1.0/TC 

)1=3.0*POLE 

)2«=3.0*POLE**2 

)3=POLE**3 

L=BB1 

L-BB2 

L=-ADl-(AAll+AA22) 

2=(BBl*AA22-BB2*AAl2) 

2=(BB2*AAll-BBl*AA21 ) 

L=AD3/(BB2*AAll-BBl*AA21) 

2=AD2-(AAll*AA22-AAl2*AA21 )+BB2*Kl 

2-(Cl*B2-C2*Bl)/(Al*B2-A2*Bl) 

3=(C2*A1-C1*A2)/(A1*B2-A2*B1) 

IME-100.0 

r=o.i 

-TIME/DT 
Y=TL/DT 

SI=  0.0 

0.0 

0.0 

=  0.5 

=  0.0 

COUNT=0 

CON=Y 

0  10,  1=1, N 
T=  T+DT 

PSIDOT=R 

VDOT   =AAll*V  +  AA12*R  +  BBl*DEL 
RDOT   =AA21*V  +  AA22*R  +  BB2*DEL 
XDOT   =COS(PSI )-V*SIN(PSI ) 
YDOT   =SIN(PSI )+V*COS(PSI ) 

PSI=PSI+(PSIDOT*DT) 
V=V+(VDOT*DT) 
R=R+(RDOT*DT) 
X=X+(XDOT*DT) 
Y=Y+( YDOT*DT) 

ICOUNT=ICOUNT+l 

IF  ( ICOUNT.NE. IY)  GO  TO  11 

YCON=Y 

ICOUNT=0 

PSI C=-ATAN ( YCON/D ) 
DEL=Kl*(PSI-PSIC)+K2*V+K3*R 
DEL=D0*TANH( DEL/DO ) 

WRITE(20,*)  T,Y 
WRITE(30,*)  T,Y 
CONTINUE 
END 
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APPENDIX  C 


PROGRAM  THESNEWT.FOR  (TWO  POINT  TIME  DELAY  ANALYSIS) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , L , NR , NV, NDRS , NDRB , I Z , MASS , 
&  NRDOT,NVDOT 
DIMENSION  A (4, 4) , FVl ( 4 ) , IVl ( 4 ) , ZZZ ( 4 , 4 ) , WR( 4 ) , WI ( 4 ) 
DIMENSION  DIST(10,10) 
OPEN  (11, FILE=' TCN8.DAT' , STATUS- 'NEW ) 


Vehicle 

i    Parameters: 

U= 

5.0 

RAT 10= 

-    0.0 

WEIGH! 

'=435.0 

IZ 

=  45.0 

L 

=  7.3 

RHO 

=  1.94 

G 

=  32.2 

XG 

=0.0104 

MASS 

=WEIGHT/G 

YRDOT  =-0.00178*0. 5*RHO*L**4 

YVDOT  =-0.03430*0. 5*RHO*L**3 

YR  =+0.01187*0. 5*RHO*L**3 

YV  =-0. 03896*0. 5*RHO*L**2 

YDRS  =+0.02345*0. 5*RHO*L**2 

YDRB  =+0.02345*0. 5*RHO*L**2 

NRDOT  =-0.00047*0. 5*RHO*L**5 

NVDOT  =-0.00178*0. 5*RHO*L**4 

NR  =-0.01022*0. 5*RHO*L**4 

NV  =-0.00769*0. 5*RHO*L**3 

NDRS  =-0.3  37*0.02  3  4  5*0.5*RHO*L**3 

NDRB  =+0.28  3*0.02  3  4  5*0.5*RHO*L**3 

DH   =( IZ-NRDOT)* (MASS-YVDOT )- 

&      ( MASS*XG-YRDOT ) * ( MASS *XG-NVDOT ) 
AA11=( ( I Z-NRDOT)*YV-( MASS *XG- YRDOT )*NV)/DH 
AA12=( ( IZ-NRDOT) *(-MASS+YR )- 

&       (MASS*XG-YRDOT)*(-MASS*XG+NR) )/DH 
AA21= ( ( MASS-YVDOT ) *NV- ( MASS*XG-NVDOT ) * YV )/DH 
AA22= ( ( MASS-YVDOT ) * ( -MASS*XG+NR ) - 

&       (MASS*XG-NVDOT)*(-MASS+YR) )/DH 
BB11=( ( IZ-NRDOT) * YDRS- ( MASS *XG-YRDOT ) *NDRS )/DH 
BB12=( ( IZ-NRDOT) * YDRB- ( MASS*XG-YRDOT ) *NDRB)/DH 
BB21=( (MASS-YVDOT) *NDRS-(MASS*XG-NVDOT) *YDRS )/DH 
BB22= ( ( MASS-YVDOT ) *NDRB- ( MASS*XG-NVDOT ) * YDRB ) /DH 

BB1=BB11+RATI0*BB12 
BB2=BB21+RATIO*BB22 

WRITE  (*,1005) 
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READ   ( * , * )      TL 

RITE  (*,1006) 
EAD  ( * , *  )       TC 

PS=1.D-10 
TMAX-100000 

L-0.01 

0  4  1=1,100 

LD  -  TL  *  L/U 

DO  1  J=l,100 
TCD»TC*L/U 
POLEC=1.0/TCD 
ADl=3.0*POLEC 
AD2=3.0*POLEC**2 
AD3-POLEC**3 
Al=BBl*U*U 
Bl=BB2*U*U 

Cl=-ADl-(AAll+AA2  2 )*U 
A2*=(BB1*AA22-BB2*AA12)*U**3 
B2=(BB2*AA11-BB1*AA21 )*U**3 
Kl=AD3/( (BB2*AA11-BB1*AA21 )*U**3) 
C2=AD2-(AA11*AA2  2-AA12*AA21 ) *U* *2+BB2*U*U*Kl 
K2=(C1*B2-C2*B1)/(A1*B2-A2*B1 ) 
K3=(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 

Al=  (AAll*BB2-AA21*BBl )*U*U*U*Kl 

A2=  (AA11*AA22-AA12*AA21 ) *U*U+ ( AAl 1 *BB2-AA21 *BBl )*U*U*U*K3- 

(AA22*BBl-AAl2*BB2)*U*U*U*K2-BB2*U*U*Kl 
A3=-(AA11+AA22)*U-(BB1*K2+BB2*K3)*U*U 
B0=  (AA11*BB2-AA21*BB1 )*U*U*U*U*K1 
B1=-(AA12*BB2-AA2  2*BB1 ) *U*U*U*Kl-BB2*U*U*U*Kl 
B2=-BB1*U*U*K1 

AQ=B2*A3-Bl 

BQ=B1*A2-B2*A1-B0*A3 

CQ=B0*A1 

DET=BQ**2-4 .D0*AQ*CQ 

IF  (DET.LT.0.D0 )  GO  TO  4 

ZTl=(-BQ+DSQRT(DET) )/(2.0D0*AQ) 

ZT2=(-BQ-DSQRT(DET) )/(2.0D0*AQ) 

IF  (ZT1.GE.0.D0)  OM=DSQRT(ZTl ) 

IF  (ZT2.GE.0.D0)  OM=DSQRT ( ZT2 ) 

ALPHA1=AQ 

ALPHA2=BQ 

ALPHA3=CQ 

BETA1=-B2 

BETA2=B0+B2*A2-B1*A3 

BETA3=B1*A1-B0*A2 

OMOLD=OM 

BETA=BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD 

ALPHA=ALPHAl*OMOLD**4+ALPHA2*OMOLD**2+ALPHA3 

BEP=5.D0*BETAl*OMOLD**4+3.D0*BETA2*OMOLD**2+BETA3 

ALP=4.D0*ALPHAl*OMOLD**3+2.D0*ALPHA2*OMOLD 

Fl=BETA*DSIN(2.D0*OMOLD*TLD)-2.D0*BETA*DSIN(OMOLD*TLD) 

F2=ALPHA*DCOS ( 2 . D0*OMOLD*TLD ) -2 . D0*ALPHA*DCOS ( OMOLD*TLD ) 
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FPl=BEP*DSIN( 2.D0*OMOLD*TLD)+BETA*2.D0*TLD*DCOS(2.D0*OMOLD*TLD) 

FP2=-2.D0*BEP*DSIN(OMOLD*TLD)-2.D0*BETA*TLD*DCOS(OMOLD*TLD) 

FP3=ALP*DCOS(2.D0*OMOLD*TLD)-ALPHA*2.D0*TLD*DSIN(2.D0*OMOLD*TLD) 

FP4=-2.D0*ALP*DCOS(OMOLD*TLD)+2.D0*TLD*ALPHA*DSIN(OMOLD*TLD) 

FPRIME=  FP1+FP2+FP3+FP4 

IF  ( FPRIME.EQ.O.DO)  STOP  1112 

IF  ( F.EQ.O.DO)  GO  TO  2 

OMNEW=OMOLD-F/FPRIME 

OMDI F=DABS ( OMNEW-OMOLD ) 

IF  (OMDIF.LE.EPS)  GO  TO  2 

OMOLD=OMNEW 

IT=IT+1 

IF  ( IT.GT.ITMAX)  STOP  1111 

GO  TO  3 

2    OM=OMNEW 

XDNUMl=DSQRT(  ( B0-B2*OM* *2 ) *  *  2+ ( Bl *OM ) *  *  2 ) 
XDNUM2=DSQRT (  ( 2 . D0*DCOS ( OM*TLD ) -DCOS ( 2 . D0*OM*TLD )  ) *  *2  + 
&   (DSIN(2.D0*OM*TLD)-2.D0*DSIN(OM*TLD) )**2) 
XDNUM=XDNUMl *XDNUM2 

XDDEN=OM*DSQRT(  ( A1-A3 *OM* *2 ) * *2+ ( A2 *OM-OM* *  3 )**2 ) 
XD=XDNUM/XDDEN 
X=XD/L 
:         DIST(I,J)=X 

WRITE(11,*)  X,TL 
:       TC=TC  +  0.015 
:     1  CONTINUE 

TL=TL  +0.02 
4  CONTINUE 
:       WRITE(11,20)  ( (DIST(I, J) , J-1,10) ,1=1,10) 
20    FORMAT(10(2X,F4.2)  ) 

1005  FORMAT  ( '  ENTER  TIME  LAG' ) 

1006  FORMAT  ( '  ENTER  TIME  CONSTANT' ) 
END 
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